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QM/MM Methods for Crystalline Defects 
Part 2 : Consistent Energy and Force-Mixing 
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Abstract 

QM/MM hybrid methods employ accurate quantum (QM) models only in regions of interest 
(defects) and switch to computationally cheaper interatomic potential (MM) models to describe 
the crystalline bulk. 

We develop two QM/MM hybrid methods for crystalline defect simulations, an energy-based 
and a force-based formulation, employing a tight binding QM model. Both methods build on 
two principles: (i) locality of the QM model; and (ii) constructing the MM model as an explicit 
and controllable approximation of the QM model. This approach enables us to establish explicit 
convergence rates in terms of the size of QM region. 


1 Introduction 

Algorithms for concurrently coupling quantum mechanics and classical molecular mechanics (QM/MM) 
are widely used to perform simulations of large systems in materials science and biochemistry BE 
fTdl [T71 [22| l32l 133] . A QM model is necessary for accurate treatments of bond breaking/formation, 
charge transfer, electron excitation and so on. However, the applications of QM is limited to systems 
with hundreds of atoms due to the significant computational cost. By contrast, MM methods based 
on empirical inter-atomic potentials are able to treat millions of atoms or more but reduced accuracy 
(more precisely, they are not transferable). QM/MM coupling methods promise (near-)QM accuracy 
at (near-)MM computational cost for large-scale atomistic simulations in materials science. 

In QM/MM simulations the computational domain is partitioned into two regions. The region of 
primary interest, described by a QM model, is embedded in an environment (e.g., bulk crystal) which 
is described by an MM model. The coupling between these two regions is the key challenge in the 
construction of accurate and efficient QM/MM methods. 

A natural question is the accuracy of QM/MM models as a function of QM region size. The 
number of atoms in the QM region is a discretisation parameter and the observables of interest should 
converge to the desired accuracy with respect to this parameter. Despite the growing number of 
QM/MM methods and their applications, few of the publications have included quantitative tests of 
the accuracy of the method and its convergence with respect to possible parameters. To our best 
knowledge, there is no theoretical analysis for QM/MM methods in the literature. 

The purpose of this paper is to initiate a numerical analysis of QM/MM methods. We develop two 
new QM/MM methods for crystalline defect simulations for which we can prove rigorous a priori error 
estimates. We use the tight binding (TB) model (a minimalist QM method) as the QM model and, 
for the MM region, construct an interatomic potential (or, forces) through an explicit approximation 
of the TB model, which is reminiscent of the force matching technique m- This approach enables us 
to establish explicit convergence rates in terms of the size of the QM region. 
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try CV47AL, UK. This work was supported by ERC Starting Grant 335120. CO’s work was also supported by the 
Leverhulme Trust through a Philip Leverhulme Prize and by EPSRC Standard Grant EP/J021377/1. 
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Our analysis is based on two key preliminaries: the “strong locality” of the (finite temperature) 
tight binding model [6] and the decay estimates of equilibria in lattices with defects [Si- 

Outline. In Section]^ we review the existing QM/MM methodology for material systems. In Section 
3, we review the tight binding model for crystalline defects which we use as the QM model. In Section 
4 and we construct QM/MM coupling schemes with rigorous error estimates based, respectively, 
on energy-mixing and force-mixing principles. Finally, we summarise our findings and make some 
concluding remarks concerning practical aspects which we will pursue in forthcoming work. 


Notation. We will use the symbol •) to denote an abstract duality pairing between a Banach space 
and its dual. The symbol | • | normally denotes the Euclidean or Frobenius norm, while || • || denotes an 
operator norm. For the sake of brevity of notation, we will denote by A\a, and {b — a \ b £ A} 

hy A — a. 

For a differentiable function /, V/ denotes the Jacobi matrix. For E £ C‘^{X), the first and second 
variations are denoted by {6E{u),v) and {5‘^E{u)w,v) for u,v,w £ X. For higher variations, we will 
use the notation 5^E{uo)[ui, • • • , Uk], and d^E{uo)[u^’‘] for abbreviation when ui = ■ ■ ■ = Uk = u. 

For j G N, gf G and V £ , we define the notation 




d^V(g) 


dg, 


pi’ 


,dg, 


for p= (pi,--- ,pj) £ AE 


Pj 


The symbol C denotes generic positive constant that may change from one line of an estimate to 
the next. When estimating rates of decay or convergence, C will always remain independent of the 
system size, of lattice position or of test functions. The dependencies of C will normally be clear from 
the context or stated explicitly. 

Acknowledgements. We thank Noam Bernstein, Gabor Csanyi and James Kermode for their helpful 
discussions. The work presented here is related to ongoing joint work. 

2 QM/MM coupling methods 

The many different variants to couple QM and MM models can be broadly divided into two categories: 
energy-mixing and force-mixing. Energy-mixing methods define an energy functional that involves 
mixture of QM and MM energies, and the solution is obtained by minimizing the energy functional. 
By contrast, force-mixing methods define a system of force balance equations, where the forces involve 
contributions from QM and MM models and are non-conservative (i.e., they are not compatible with 
any energy functional). In the present section we review both classes of QM/MM schemes for materials 
systems. 


2.1 Energy-mixing 

The system under consideration is partitioned into QM and MM regions (see Figure [2T] for two exam¬ 
ples). Let and denote the respective atomic configurations in these two regions. Depending 
on the construction of the hybrid total energy E, energy-based methods can mainly be divided into 
two categories: 

(1) In the subtractive approach, e.g. the OiV/OMmethod [THl|22l |29] and its derivatives. 




( 2 . 1 ) 


that is, an MM energy for the entire system is corrected by the difference between the QM and MM 
energies of the QM region. 
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Figure 2.1: Partition of QM and MM regions for an edge dislocation in a 2D triangular lattice and 
crack in the 2D hexagonal lattice (cartoons). The dislocation core and a small neighbourhood belong 
to the QM region (red / dark), while the bulk crystal behaves purely elastically and can therefore be 
well described by an empirical interatomic potential (blue / light). 


(2) In the additive approach, e.g. ChemShell [20], DL-FIND [T6|, MAAD |3] and QUASI |27jl the 
QM energy of the QM region and MM energy of the MM region are connected via an interaction 
energy which may depend on an interface that involves parts of both regions. 


£:(yQMUMM) _ 


,MM\ 


+ E 


linteraction/ QM MM 




( 2 . 2 ) 


The advantage of energy-based methods is that they are naturally energy conserving. Unfortu¬ 
nately, the spurious interface effects acting between the QM and MM regions can be significant. To 
alleviate such effects, the QM and MM regions are either “passivated” [211 EH] or “buffered” [2ll|22|. 
In the first approach, the energies _EQ^(yQ^) in (2.1) and (2.2) are the energies of the passivated 


cluster of the QM region, in which a number of additional atoms that have no counterparts in the real 
system (for example, hydrogen atoms) are added to the QM region to terminate the broken bonds. 
The second approach handles the boundary by defining buffer layers surrounding the QM and MM 
regions, so that each atom can see a full complement of surrounding atoms. 

The second approach seems to be preferred in solid state systems since the elimination of the 
boundary effects for passivated atoms will not be perfect and may indeed be severe [7] . The simplest 
example is for a perfect bulk system, where the true force on all atoms are zero. However, the 
passivated cluster force computed with QM and MM will in general be non-zero on the passivation 
atoms and nearby atoms [Hj. This is reminiscent of the ghost forces which are a well-understood 
concept in atomistic/continuum multi-scale methods |I9j . 


2.1.1 An idealized hybrid model. Beside the widely used subtractive and additive approaches, 
there is a third type of energy-based formulations: the local energy approach, which mixes local energies 
computed by QM and MM methods in their respective regions. 


Eiy 


QMUMM^ _ 




(2.3) 


teQM 


teMM 


with E£ denoting the local energy associated to the Uth atomic site. Even though the expression (2.3) 


seems intuitive, this variant is not commonly used in QM/MM coupling schemes. Indeed, we are only 
aware of a brief reference to this approach in [3| . The reason is that it was unclear how to decompose 


E into local contributions E^^ that match with a classical interatomic potential site energy E^ 


MM 
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In the previous work of this series [H] , we studied the tight binding site energy introduced in mm 
and justified its “strong locality” rigorously. This is important and useful in QM/MM coupling scheme 


based on (2.3) since: (1) when using classical potentials the total energy is almost always written as a 
sum over atoms E = E^, therefore, we are able to establish the bridge between the electronic and 

classical worlds; (2) rather than using “black-box” MM potentials, we can construct MM site energies 
based on the approximations of QM site energies for good coupling of the different models. 

It is pointed out in [3] that matching the force-constant/dynamic matrix (i.e. the first derivatives 
of the force or the second order derivatives of the energy with respect to atomic positions) would 
guarantee a perfect match between the MM and QM forces for arbitrary infinitesimal displacements 
from equilibrium. In case of an energy-based method, only such a strict matching criterion can 
guarantee that spurious forces are eliminated for near equilibrium configurations. The errors resulting 
from mismatching the force-constant/dynamic matrix are analogous to so-called “ghost forces” in 
atomistic/continuum hybrid schemes. 

Based on these observations, we construct an idealized MM site energy by taking the second order 
expansion of 


■= (yo) + {SEl {yo),y - yo) + -{6'^E^^{yo){y - yo),y - yo) 


(2.4) 


where y = -g ^ predicted (near-)equilibrium configuration, typically the far-held 


crystalline environment or an explicit linearlised elasticity solution. The QM/MM total energy (2.3) 
is then given by 


E{y) = 


E 

teQM 


piQM 


iv) + E 


(2.5) 


£eMM 


(2.5) gives rise to a simple QM/MM coupling scheme, in which the MM potential is constructed such 
that it matches the QM model. Matching of the MM and QM models for higher order information can 
also become important, e.g., for slowly decaying elastic helds (dislocations, cracks) or due to increased 
temperature which may cause huctuations to displacements beyond the quadratic regime [3]. We will 
therefore discuss arbitrary order expansions in this paper. 


2.2 Force-mixing 

With the partition of QM and MM regions, the force-based methods combine QM forces for atoms 
in the QM region with MM forces in the MM region. The simplest variant, brutal force mixing m , is 
defined by 


i?^(yQMUMM) 


^QM(yQMUMM) if £ G QM 


( 2 . 6 ) 


Typically, the QM and MM forces are computed by carving a cluster buffered by a layer of atoms 
defined by a distance cutoff, and the forces on all the buffer atoms are discarded [7j . In other variants, 
a transition region is introduced where forces are smoothly blended [7]. Examples of force-mixing 
QM/MM methods are DCET and LOTF 

The main advantage of force-mixing is that there are no spurious interface forces as in energy¬ 
mixing schemes. However, this comes at the cost of a non-conservative force held. Moreover, if the 
QM and MM forces are directly used (without modihcation) for molecular dynamics simulations, then 
the dynamics will not conserve momentum [3]. 
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2.2.1 An idealized model. Similar to the discussions for energy-mixing methods, we will construct 
MM forces by an expansion of QM forces so that the force-constant/dynamic matrix can be matched, 
e.g., with a first order expansion, 

F^^iy) := Ff^iy) := F^^{yo) + {6F^^iyo),y - yo), (2.7) 

where y = yQ^uMM jg suitable predictor of the equilibrium configuration. For the same 

reasons as in the energy-mixing approach, we will also consider higher-order expansions of the forces. 

Remark 2.1. Our construction of the MM site energies and of the MM forces is reminiscent of the 
classical idea of force matching. This is usually applied to the construction of interatomic potentials 
m and has more recently been applied in a coupling context, e.g., in BW- A key difference in our 
present work, in the energy-based variant, is that we match the site energies rather than that total 
energies (and forces). 


3 Tight binding model for crystalline defects 

A finite or countable index-set A C M”* is called a reference configuration. A deformed configuration 
is described by a map y : A ^ with m,d ^ {2, 3} denoting the space dimensions. (Allowing m d 
allows us to define 2D models of straight dislocations.) 

We say that the map y is a proper configuration if the atoms do not accumulate: 

L. 3 m > 0 such that \y{A) — y{k)\ > m|^ — k\ V £, /c G A. 

Throughout, we let Vm C (M'^)^ denote the subset of all y G (M'^)^ satisfying L. If we need to emphasize 
the domain A, then we will write Vm(A). 


3.1 The tight binding model and its site energy 

The tight binding model is a minimalist QM type model, which enables the investigation and prediction 
of properties of molecules and materials. For simplicity of presentation, we consider a ‘two-centre’ 
tight binding model [HIES] with a single orbital per atom and the identity overlap matrix. All results 
can be extended directly to general non-self-consistent tight binding models, as described in |6|. 

For a finite system with reference conhguration D, ffkl = N, the ‘two-centre’ tight binding model 
is formulated in terms of a discrete Hamiltonian, with the matrix elements 



if 7 = A: 
if i^k, 


(3.1) 


where Rc is a cut-off radius, /ions G C'’’([0, oo)) is the on-site term, with g G C'"([0,oo)), ^ = 0 in 
[Rc,oo), and /ihop G C'’’([0, oo)) is the hopping term with /ihop(?’) = 0 Vr G [Rc,oo). 

Our results can be generalised to the more general TB model presented in |6|, but for the sake of 
simplicity of notation, we restrict ourselves to (3.1), which still includes most non-self-consistent TB 
models in the literature. 

Note that /ions and /ihop are independent of i and k, which indicates that all atoms of the system 
belong to the same species. We observe that the formulation (3.1) satisfies all the assumptions on 
Hamiltonian matrix elements in jh] Assumptions H.tb, H.loc, H.sym, H.embj. 

With the above tight binding Hamiltonion 77, we can obtain the band energy of the system 


N 

E^iv) = '^f{^s)es 

s=l 


(3.2) 
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where are the eigenvalues of T-L{y), with associated eigenvectors V's, 


/ the Fermi-Dirac function, 


'Hiy)i>s = Ssi’s s = l,2,---,N, 


f{e) = (l + ^ ^ 


(3.3) 


(3.4) 


/i a fixed chemical potential, k-Q Boltzmann’s constant, and T > 0 the temperature of the system. We 
do not consider the pairwise repulsive potential, which can be treated purely classically |6j. 

Following |T2], we can distribute the energy to each atomic site 


E^{y) = '^E^{y) with E^{y) :='^f{es)e,\['il;s]e\‘^ ■ 
e&n s 


(3.5) 


The following theorem [6l Theorem 3.1 (i)] states the existence of the limit as f A, C A, for 
some countable reference domain A. For an infinite body, A, we will denote this limit site energy by 
Ei. We will continue to denote the site energies of subsystems H C A by E^. When A is an infinite 
reference configuration and A C then we will also use the short-hand Ef := Ef^^. 

Theorem 3.1. Suppose A is countable, y G Vm(A) a deformation, and H C A a finite subset. Then, 

(i) (regularity and locality of the site energy) Ef{y) possesses jth order partial derivatives with 


1 < J < n — 1, and there exist constants Cj and rjj such that 


d^E^{y) 


< Cje~^^ 


(3.6) 


d[y{'mi)]n ■ ■■d[yimj)]i 
with mk G n and 1 < ffc < d for any 1 < k < j; 

(ii) (isometry invariance) Ef‘{y) = E^{g{y)) if g : M'’* —>■ M”* is an isometry; 

(in) (permutation invariance) E^{y) = {y o Q) for a permutation Q : A —)• A. 

(iv) (thermodynamic limit) Efiy) := lim/j_^oo exists and satisfies (i), (ii), (in) with H = A 

For a finite subset fl C A, we define the (negative) force 


E^\y) :=VE^fiy), 


m 


M dE^(y) 

component notation, [E^ {y)\. = 


1 < i < d. 


Using (3.5), we have 


[T{!/)], = E 


fcen 


dE^jy) 
d[yii)h ’ 


(3.7) 


(3.8) 


which, together with Theorem |3.1t yields the thermodynamic limit of the force, as well as its regularity, 
locality, and isometry/permutation invariance. 

Theorem 3.2. Suppose the assumptions of Theorem\3.1\ are satisfied, then 


(i) (regularity and locality of the force) Ep{y) possesses jth order partial derivatives with 1 < j < 
n — 2, and there exist constants Cj and gj such that 


d^FPiy) 


5 [ 2 /(" 1 i)]u • ■■d[yimj)]i 
with mk G U and 1 < ffc < d for any 1 < k < j; 


< Cje~^^ 


(3.9) 
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(ii) (isometry invariance) Fp(^Qy + c) = Q^Fp{y) for all Q G SO{d),c G 
(Hi) (permutation invariance) Fp{y) = Fg_il^^\y o Q) for a permutation G : A ^ A; 

(iv) (thermodynamic limit) -F)(y) := limij^oo Ff^^^py) exists and satisfies (i), (ii), (Hi) with Gl = A. 


3.2 Crystalline defects 

3.2.1 Energy space. Let A C be an infinite reference configuration satisfying A \ = 

(AZ”*) \ -Brdef where iinEF > 0, A G SL(m) and A n -Bhdef finite. For analytical purposes, we 
assume that there is a regular partition 7a of into triangles if m = 2 and tetrahedra if m = 3, 


whose nodes are the reference sites A (see Appendix A for interpolations of lattice functions on this 
background mesh). 

We can decompose the deformation 


y(^) = Uoi^) + ui^) = Pi + uo{£) + u{i) y i € A, 


(3.10) 


IS a 


where uq : A —>■ is a predictor prescribing the far-field boundary condition, u : A —)> 

corrector, and P G denotes a macroscopically applied deformation. 

If G A and i + p € A, then we define the finite difference Dpu{£) := u{£ + p) — u{i). If 7^ C A — 
then we define D'jiu{l) := {Dpu{£))p^'ii, and Du{£) := For a stencil Du{i) and 7 > 0 we 

define the (semi-)norms 

E \\Du\\, 2 := (Y,\Dnie)\U . 

'^peA-£ 2 £eA 2 


An immediate consequence of (A.l) is that all (semi-)norms 


!, 7 > 0 , are equivalent. 


We can now define the natural function space of finite-energy displacements, 

;= ; A ^ \\Du\\i 2 < 00}. 


3.2.2 Site energy. Let Ei denote the site energies we defined in Theorem 
are translation invariant we can define lA • (M™)^”^ —?• M by 


3.1 


iv). Because they 


VpDu) := EpPxo + u) with xq : A 


and xo(£) = i y i £ A. 


(3.11) 


For a displacement u satisfying y^ + u ^ Vm (A), we can define (formally, for now) the energy-difference 
functional 

- E + '^) - ^^( 2 ^ 0 )) = E {y^iDuo{£) + Du{£)) - VPDuo{£))) . (3.12) 


£eA 


teA 


For both point defects and dislocations, we can construct predictors 7/0 (see ^3.2.3 and ^3.2.4) 


such that (IT(0) G . We prove in [5] (see also [9l Lemma 2.1]) that, under this condition, £ is 

well-defined on the space Admo and in fact £ G C"“^(Admo), where 

Adm^ := {u G \yQ{£) + u{£) — yQ{m) — u{m)\ > m\i — m\ V£,mGA}. 


In ^3.2.3 and ^3.2.4 we show how the crucial condition (5T(0) G is obtained for, respectively, 


point defects and dislocations. In § 3.2.5 we then present a unified description for which we then 
rigorously state the properties of £ and the associated variational problem. 
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3.2.3 Point defects. We make the following standing assumptions for point defects: 

P. m = d e {2, 3}; 3 Rbef > 0, ^ G SL(m) such that and A n -B^def i® 

finite; P = Id; uq = 0; yo{i) = I- 


3.2.4 Dislocations. The following derivation is not essential to our analysis of QM/MM schemes, 
and can indeed be found in [9], however, for the sake of completeness, we still present enough detail 
to justify the unified formulation in § 3.2.5 

We consider a model for straight dislocation lines obtained by projecting a 3D crystal. For a 3D 
lattice Bl? with dislocation direction parallel to 63 and Burgers vector b = ( 61 , 0 , 63 ), we consider 
displacements W : BT? —)> that are periodic in the direction of the dislocation direction 63 . Thus, 
we choose a projected reference lattice A := Al? = {(£i,£ 2 ) \ ^ = (-^ 1 ,^ 2 , ^ 3 ) £ which is again a 

Bravais lattice. This projection gives rise to a projected 2D site energy with the additional invariance 


Ei{y) = Ei{y + 263 ) Vz : A b 3 Z 


(3.13) 


Let X G be the position of the dislocation core and T := {x G | X 2 = X 2 , xi > xi} be the 
“branch cut”, with x chosen such that TnA = 0. Following [9], we define the far-held predictor uq by 


where G C“(M2 \ P; 


uo(x) := n'“(e(x)), (3.14) 

is the continuum linear elasticity solution (see [9] for the details) and 


^(x) = X - bi2—?? 


X — X 


arg(x — x). 


(3.15) 


with 


with arg(x) denoting the angle in (0, 27r) between x and bi 2 = (bi, b 2 ) = (bi, 0), and r/ G 
r/ = 0 in (— 00 , 0 ], r/ = 1 in [ 1 , 00 ) removes the singularity. 

The predictor yo = Px + uq is constructed in such a way that yo jumps across T, which encodes 
the presence of the dislocation. But there is an ambiguity in this dehnition in that we could have 
equally placed the jump into the left half-plane {xi < xi}. The role of ^ in the dehnition of uq is that 
applying a plastic slip across the plane {x 2 = X 2 } via the dehnition 


/(x) := 


y(£), 


h > X2, 


y(£ - bi2) - b3e3, h < X2 


achieves exactly this transfer: it leaves the (3D) conhguration invariant, while generating a new 
predictor y^ G (^“(Dr) where Dp = {xi > xi-|-x-|-bi}. Since the map y ^ y^ represents a relabelling 


of the atom indices and an integer shift in the out-of-plane direction, we can apply (3.13) and Theorem 
3.1| (hi) to obtain 

Ei{y) = Es*e{y^), (3.16) 

where S is the ^^-orthogonal operator with inverse S* = S~^ dehned by 

Suit) := I “!(>' , , and S'u(l) := | “<(>■ , , 

^ ^ \ u{£ - bi 2 ), £2 < X2 ^ \ tt(£ -I- bi 2 ), £2 < X2. 


We can translate (3.16) to a statement about uq and Let Sqw{x) = w(x),X 2 > X 2 and 
Sow{x) = w{x — bi 2 ) — b,X 2 < X 2 , then we obtain that y^ = Px P Squq and Squq G (^“(Dr) and 


So{uo + u) = SqUq + Su. The permutation invariance (3.16) can now be rewritten as an invariance of 
Vi under the slip Sq: 


V{D{uo + u){£)) = V(e(£) + Du(£)) Vn G Admo, £ ^ A 


(3.17) 
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where 


e{t) := (ep(£))pgA_£ with ep{t) := | 


S*DpSou^{i), 

DpUo{i), otherwise, 


and 


Du{£) := {Dpu{£))p(zA-e with Dpu{£) := 


S*DpSu{i), i£nr, 
Dpu{i), 


otherwise. 


(3.18) 

(3.19) 


The following lemma, proven in [5], is a straightforward extension of [HJ Lemma 3.1]. 


Lemma 3.1. If the predictor uq is defined by (3.14) and e{t) is given by (3.18), then there exists a 
constant C such that 


|eo-(^)| < C|cj| • 1^1 ^ and \Dpea{I)\'Y < C\p\ ■ |(t| • |.^| 

We now summarise our standing assumptions for dislocations: 

D. m = 2, d = 3; A = AZ”*; P{£i,i 2 ) = ^ 3 ); uq is defined by ( 3.14[ ); yo{i) = £ + uo(^)- 


(3.20) 


Remark 3.1. One can treat anti-plane models of pure screw dislocations by admitting displacements 
of the form uq = (0, 0 ,?xo, 3 ) o,nd u = (0, 0 , 7 / 3 ). Similarly, one can treat the in-plane models of pure 
edge dislocations by admitting displacements of the form uq = (uo,i, uo, 2 ,0) and u = {ui,U 2 , 0 ) fSi- 
For anti-plane models the atoms do not accumulate and the condition L can be ignored. 


3.2.5 Unified formulation. In order to consider the point defect and dislocation cases within a 
unified notation we introduce the following notation. Let 



if P 
if D 


e(£) := 


0 if P 

e{i) ifD 


and Du{l) := 


Du{I) ifP 
Du{e) ifD 


Using the assumption uq = 0 in P for point defects and the slip invariance condition (3.17) for 
dislocations, we can rewrite the energy difference functional (|3.12[) as 


£{u) = Y, + Du(£)) - W(e(^))) 

7eA 


(3.21) 


for both point defects and dislocations. The following result is proven in [5|, extending an analogous 
result in [9] which is restricted to finite-range site energies. 

Lemma 3.2. Suppose that P or D is satisfied, then E is well-defined on n Admo, where 

:= {tt : A —>■ I 3R > 0 s.t. u = const in A \ 

and continuous with respect to the -topology. Therefore, there exists a unique continuous extension 
to which belongs to C'"“^(l^^’^). 

Having a well-defined energy-difference functional, the equilibrium state can be determined by 
solving the variational problem 


u G argmin |T(m), M G Admo}, (3.22) 

where “arg min” is understood in the sense of local minimality. 

The next result is an extension of P Theorem 2.3 and 3.5], which gives the decay estimates for 
the equilibrium state for point defects and dislocations (see [3] for the proof). 
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Theorem 3.3. Suppose that either P or D is satisfied. If u £ Admo is a strongly stable solution to 

3 c > 0 s.t. (u)r, r) > c||Dr ||^2 Vr G 


(3.22), that is, 


then there exists a constant C > 0 and Uoo G such that 

<c{i + \i\)-'^\og\2 + \i\), 

\u{i)-Uoo\<C{l + \l\Y-^\og\2 + \I\), 
where t = 0 for case P and t = 1 for case D. 


(3.23) 

(3.24) 

(3.25) 


Remark 3.2. It can he immediately seen that \Du{i)\.y satisfies the same estimate as (3.24). 


4 Energy-mixing 

4.1 Formulation of QM/MM energy mixing 

Following the outline in § 2.1.1[ we construct approximations to the tight binding site energy V£{g) 


y^MM(gj) g g (]^™-)^ t i^y Taylor’s expansion, and approximate the energy difference functional by 
S{u)^ E {Ve{e{£) + Du{i)) -Ve{B{£))) + E {e{i) + Du{e)) - {e{£))). (4.1) 

£gAQM 


£eAMM 


Since minimising (4.1) over u G Admo is an infinite-dimensional problem, we will also approximate 


the space of trial functions. 

4.1.1 Decomposition of A. We decompose the reference configuration A into three disjoint sets, 
A = A^'^ U A'^^ U A^^, where A^^ denotes the QM region, A^^ the MM region and A^^ the far-field 
where atom positions will be frozen to those given by the far-held predictor. Moreover, we dehne a 
buffer region A^^^ C A'^^ surrounding A*^^ such that all atoms in ABUf 

are involved in the tight 

binding calculation when evaluating the site energies in A^'^. 

For simplicity, we use balls centred at the point defect or dislocation core to decompose A, and 
use parameters Rqm, Rmm and Rbuf to represent the respective radii, that is, 

= SitQM n A, A^^^ = n A \ aQ^^, aI'uf ^ (^Br^^+r^^^ n A) \ aQ^, 


and A^^ = A \ (A^^ U A'^^). See Figure 4.2 for a visualisation of this decomposition 


4.1.2 Buffered QM model and site energies. The site energies in the exact model have inhnite 
range, henve we truncate them to obtain a computable approximation. To that end, we dehne 




BUF, 


g ■= 








[9h 


I G A^^^^ U A^^ 


(4.2) 


where Vp{Du{l)) := E^{Pxq + u). 

We assume throughout that Rqm > Rdef + Rbuf- In this case, Theorem |3.1| (ii) (iii) and the 
assumptions on A in P and D imply that the truncated site potential (4.2) is independent of I in 
AMM y A^^. That is, there exists : (M”*)^ —>■ M such that 

= Vf^^{Dnu{e)) with IZ = Br^^^ n (A \ 0), ^ i£ A^^ U A^^. (4.3) 
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Figure 4.2: Decomposition of a crystal lattice with defect (edge dislocation) into QM, MM, buffer 


and far-held regions, according to § 4.1.1 


Remark 4.1. IFe have used the buffer radius parameter Rruf for both the buffer surrounding the 
QM region and for the buffer used in the approximate site potential . Although we could choose 
two separate parameters, they affect the error in similar ways, hence for simplicity of notation we use 
only one parameter. 


4.1.3 QM/MM coupling. The homogeneity (4.3) allows us to construct the MM site potential 
by A:-th order Taylor expansion of about the far-held lattice state. 


1 


_ TkVf^^{g) := [g®^] with k>2. 

1=1 


(4.4) 


With the dehnitions (4.2) and (4.4) we can now specify the QM/MM energy-mixing scheme 




G argmin|T^(u) | u G Adm^j, 


(4.5) 


with the QM/MM hybrid energy difference functional 

aBUF 


E [vP''^{e{i) + Du{l))-V^ 


BUF/ 




^eA^MuAi^r’ 

and admissible set 

AdmJJ := AdiUm n where 


^ ('yMM ^ _ yMM 


(4.6) 




'W^ ;= <! M G 


u = 0 in A^^l 


(4.7) 


Using same arguments as those in lai, we have that G C" ^(Adm^). Note, in particular, 
that the sum over A^'^ U A^^ in the dehnition of is in fact hnite, since has a hnite range of 

interaction. 
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4.2 Error estimates 


The QM/MM energy mixing scheme (4.5) satisfies the following approximation error estimate. The 
main steps of the proof are presented below, but some technical details are given in the appendix. 

Theorem 4.1. Suppose that either assumption P or D is satisfied and that u is a strongly stable 
solution of (3.22). 


If, in the definition of in (4.6), is the k-th order expansion in (4.4) and n > k + 2, then 

there exist constants C, k, Cg^p, such that, for Rqm sufficiently large and for 

-Rbuf > max{c§up log Rqu, Cg^ log log Rum} 


there exists a strongly stable solution of (|4.5l) satisfying 


\\Du- Du^\\i 2 < C (^Rqm + W -Rmm + e ^ 

\£{u) - £^iu^)\ < C (4?qm log* ^QM + Rmm log'* ^mm + 


—k/Jbuf 


(4.8) 

(4.9) 


where 


a = {2k — l)m/2, i/P, 
a = k — 1, */ D. 

Proof. 1. Quasi-best approximation: Following [9l Lemma 7.3], we can construct T^u G AdmQ by 

T^u{i) :=r]{i/RMM){u{i)-Uoc-aR^^), (4.10) 


where Uoo is given in Theorem 


3.3 




•- ih 


m 


^Appendix A and rj G C'^(M™') is a cut-off function satisfying r]{x) = 1 for jxj <4/6 and r]{x) = 0 for 
|x| >5/6. Then, for any 7 > 0 and for i?MM sufficiently large, we have from the decay estimates in 
Theorem 13.31 that 


|L>r“n- 


u - Du \\^2 < and 

\DT^u{e)\-^ < C(l + K|)-”"log*(2 + |£|) V^G A. 


(4.11) 

(4.12) 


Let r > 0 be such that Br{u) C Adm^ for some m > 0. We have from Theorem |3 . 3| that, for i?MM 
sufficiently large, T^u G Bj./ 2 {u) and hence Bj./ 2 {'R^I'-) C Admn,. Since £ G C^(Admo), 6£ and d‘^£ 
are Lipschitz continuous in Br{u) with uniform Lipschitz constants Li and L 2 , i.e.. 


\\5£{u) - 5£{T^u)\\ < Lfi\Du-DT^u\U 2 < and (4.13) 

\\5‘^£{u) - 5^£{T^u)\\ < L 2 \\Du - < CL2\\Du\\,2^R^B^^^/,y (4-14) 

2. Stability: Since u is strongly stable, there exists c > 0 such that {S‘^£^{u)v,v) > c||Iln||rt- For 
any v G we have 

{6'^£^{T^u)v,v) — (6‘^£{u)v,v) 

= {{6'^£^{T^u) - 6‘^£{T^u))v,v) + {{6^£{T^u) - 5^£{u))v,v) 

= ^ ^ DT^uiQ) - S'^Vi{e{£) + DT^u{i))'^Dv{£), Dn(£)' 


eeA 

+ 


E 

£eAMMuABF 


+ DT^u{i)) - + DT^u{e))'^Dv{i), Dn(^) 


+ {{6^£{T^u)-6^£{u))v,v) 

Qi + Q2 P Qs- 


(4.15) 
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Using the estimate (B.3), we have 


£eA 


Taylor’s expansion (4.4) yields 


(4.16) 


l^2l = + DTM^)) - 5V|UF(e(t’) + DTMi)))Dv{i), Dt;(£)) 

^eAMMu^FF 

<C \e{i) + DTM^t^~"\Dv{ity < + (4.17) 

feAMMuAFF ^ 


The Lipschitz continuity (4.14) implies 


IQsl < CL2 ||M||,;|a\b,„„„)P«-II,V 


(4.18) 


Using (4.12), (4.15), (4.16), (4.17), ( |4.18[ ), the decay estimates in Lemma [3T| Theorem 3.3 and the 


fact that M is a strongly stable solution, we have that for sufficiently large Rqm and iiBUF (note that 
.Rmm > Rqm), 


{6^S^{T^u)v,v)>^\\Dvf,2. (4.19) 

3. Consistency: We estimate the consistency error, for any v G by 

{SS^{T^u),v) 

= {5£^{T^u) - 6£iT'^u),v) + {5£{T^u) - 6£{u),v) 

= + DTH^(^)) - 5Ve{e{£) + DT^u{£)),Dv{i)) 

ieA 

+ (5U^^(e(£) + DT^u{i)) - ^ DT^u{i)),Dv{£)) 

ieA^Mu^FF 

+{ 6 £{T'^u) - 6 £{u),v) 

:= T 1 + T 2 + Tg. (4.20) 

The term Ti can be estimated by 


|ri|<Ce-««BUFp^|| 


with some constant k; a detailed proof of this assertion is presented in Appendix B 
To estimate T 2 , we have from (14.4^ that 


(4.21) 


feAMMuAFF 

< C \e{£) + DTM£)\’;\Dv{i)U 

^eAMMuaff 

< Clle + DT^-u||^2fc(AMMuAFF)ll-D^ll£2 (4.22) 
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Further, using (4.13) we can estimate by 


IT 3 I < CLi||F)u||£ 2 (a\b^^^^ 2 )||F>u||£ 2 . (4.23) 

Taking into accounts (4.20), (4.21), ( 4.22[ ) and ( |4.23[ ), we have 

{SS^{T^u),v) < C (e-"«BUF + ||e + DrHu||,^,,(^MMuAFP) + \\Du\\mA\Bn^^/,)) \\Dv\U2- (4-24) 


If P is satisfied, then we can obtain the estimates for point defects by substituting e(i) = 0 and 
\0u{i)\^ < C{1 + into ( |434l ): 




(4.25) 


IfD is satisfied, then we can obtain the estimates for dislocations by substituting e{£) = e{£), |e(£)|.y < 


C\e\-\ and \Duie)\^ < C (1 + |^|)-2 log (2 + |£|) into ( [4^ : 

\{6£^iT^u),v)\ < C + R^^logRuu + 


\\Dv\U- 


(4.26) 


4- Application of inverse function theorem: With the stability (4.19) and consistency (4.24), we 
can apply the inverse function theorem [23l Lemma 2 . 2 ] to obtain, for Rqm, Rbvf sufficiently large, 
the existence of a solution to (4.5), and the estimate 


\\Du^-DT^u\\p<C 


'H-i 


J} l)m/2 jj m/2 -r p 


^QM^^ + -^MM log -RmM + e if D 


(4.27) 


which together with (4.11) completes the proof of (4.8). The error estimate, together with the stability 
estimate (4.19) in particular imply that, for Rqm, Rb\jf sufficiently large, is strongly stable. 

5. Error in the energy: Next, we estimate the error in the energy difference functional. From 
£ G C'^(Admo) we have that 

\£{T^u)-£{u)\= [ {6£{{1 -s)u + sT'^u),T'^u-u) ds 

Jo 

f {5£{{1-s)u +sT'^u) - 5£{u),T'^u-u) ds < C\\DT'^u - DuW'j^ 


*7 


and from £^ G C'^(Adm^) that 

|£:^(u^) - £^{T^u)\ < C\\Du^ - DT^ 

Denoting g{£) = DT^u{i) and suppressing the argument {£) in g{£) and e{£), we have 

\£{T^u) - £^(T^u)\ 

= {^^(9 + e) - Ve{e) - V^^^ig + e) + 

i&k 

+ Y. + e) - y|U^(e) - nV^^^{g + e) + 

^sAMMuaff 


(4.28) 


(4.29) 


:= S 1 + S 2 , 

where 5i is estimated in [Appendix B|by 


|5i| < Ce"^^BUF^ 


(4.30) 

(4.31) 
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and S 2 is estimated by 


1^21 < 


^sAMMuAr’^’ 


1 


{k + l)l 




+C(|5 + e|^+2^|e|^+2 


< 


^ X/ (1517(15 + + |e|^) + 1 ^ + 


feAMMuAFF 

<c(mT^u( 


<c 


R. 


— km 
QM 


1^ • |e(£) + DT tt(^)|^||£i(AMMuAFF) + |||e(£) + DT 

if P 


k+2 


^QM log ^QM if D 


Taking ( |4.27| ), ( 4.28| ), ( 4.29| ), (4.30), (4.31) and (4.32) into accounts, we have 


£^{u^)-£{u)\ < \£^{u^)-£^{T^u)\ + \£^{T^u)-£{T^u)\ + \£{T^u)-£{u 




7 H/TnH- 




nH- 


nH- 


(4.32) 


< C 


Tj—km _|_ I p—nR^up 

' -'•'A/TM W c 


■-MM 


if P 


log Rqm + Rmm log^ Rmm + e if d 


which completes the proof of (4.9). 


(4.33) 


□ 


5 Force-mixing 

5.1 Formulation of QM/MM force mixing 


To construct a force-based QM/MM coupling scheme, we follow the idea in § 2.2.1 In the MM region 
we construct an approximation to the tight binding force ^^(y) ~ by Taylor’s expansion in 

order to ensure a good match between the QM and MM models. 


Our starting point, instead of the energy minimisation formulation (3.22), is the force-equilibrium 
formulation: 


Find u G Admo, s.t. Fiiyo -|- u) = 0 V £ G A. 


(5.1) 


where 


Fi{yo + u)= Y, V,_p^p[Duo{^ - p) + Du{t - p)) - Y Vi^p{Duo{^) + Du{e)). (5.2) 

p^l-k p&k-i 


We have from (3.12) and Theorem |3.2| (iv) that F(yo + u) = V£{u), hence any solution of (3.22) also 
solves (5.2). 

To simplify the notation in the construction of the QM/MM scheme we define 
Fp{u) := Fp{yo + u) and F^iw) := Fp{Pxo + w), 

and we remark that 

FP{uo + u)=fP{u)= Y vY,p{DMe-P) + ^<^-P))- E VPp{Duoii) + Du{i)). (5.3) 

p&i-fi p&fi-i 
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We decompose the reference configuration into A^^, A®^^, in the same way as in §|^ 

To obtain computable forces we then truncate the force of the infinite lattice, 


: = 


£eAQM, 

j-SflBUF W ^ ^ g ^MM u ^FF _ 


(5.4) 


If Rqm > .Rdef + .Rbuf, then Theorem 3.2 (ii) (iii) and the assumptions on A in P and D imply 

that the truncated force operator is independent of i in A'^^ U A^^. That is, there exists 

R, where TZ = (AZ”*) n such that 


jpBUF . ^J^rr^)7^ 




(5.5) 


y A^^ u A^^. 

We now define the MM force to be the k-ih order Taylor expansion of 

k 

:= TkRf^^iw) := T'|UF(o) + ^ J-|UF(o) with A: > 1. (5.6) 


i=i 


We remark that the zeroth-order term in the expansion vanishes since the reference lattice is an 
equilibrium. 

We have the following force-mixing QM/MM coupling model: 


Find £ Adm^, s.t. (u^) =0 y i£ A^^ U A^^ 
with the hybrid force 

e £ aQ“ 


T-f (u) = 






l£ A 


MM 


We emphasize that J-^ is not a gradient of any energy functional. For v : A 
notation (T'^(u),n) := ■ v{i) in our analysis. 


(5.7) 


(5.8) 


we will use the 


5.2 Error estimates 

Theorem 5.1. Suppose that either assumption P or D is satisfied and that u is a strongly stable 
solution of (5.1). 


Suppose that, in the definition ofT^ in (5.8), is the k-th order expansion in (5.6) andn > k+ 

3. Then, for any given MM region growth constant Cq^ > 0, there exist constants C, k, Cg^p, Cg^ 
such that, if Rqm is sufficiently large, while Rqm,Rb\jf,Rmm maintain the bounds 

log < Cqm and i2BUF > max{ c^gp log Rum , Crot log log Rum }, 

Rqm 


there exists a strongly stable solution of (5.7) satisfying 

\\Du- Du^W^ < C log 72 mm + log* 72 mm + ^ 

f a = (2A: + l)m/2, i/P, 
where < 

1 a = A;, */ D. 


(5.9) 
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Remark 5.1. In view of the bound log < Cq^ we could replace logi^MM with logi?QM ^ (5.9), 
however, we keep log Rmm to highlight the dependence of the error estimate on the growth of Rmm 
relative to Rqm- 


Proof. We will follow the same strategy as the proof of Theorem |4.1[ 

1. Quasi-best approximation: We take the approximation T^u G Adm[^ constructed in the proof 
of Theorem 4.1, so that the properties from (4.11) to (4.14) are satisfied. 

2. Stability: Let be defined by (4.6) with l/'^''^^Tdng the {k + l)-th order expansion in (4.4). 
For any v G we have 

(bT^{T^u)v,v) = (( 5 .T^(r^h) - b^E^{ff^u))v,v) + (fE^[T^u)v,v), 
where the first term is estimated in [Appendix C 


(5.10) 


as 


\{{5R^iT^u) - 5^S^{T^u))v,v)\ < (7 ||^^|| 2 ^ 


(5.11) 

with some constant k, and the second term is estimated in ^ (4.19). Therefore, we have that for 
sufficiently large Rqm, Rmm and Rbuf, 


{6T^{T^u)v,v)>-\\Dv\\j,. 

3. Consistency: We estimate the consistency error for any v G 

lQF^{T^u),v) = {bE^(T^u),v) + - bE'^{T^u),v), 

where the first term has been estimated in ^ and the second term can be written as 
{T^{T^u) - bE^{T^u),v) 

= Y. ^ [Tf{T^h) - V,T«(tHu))u(£) 

+ ^ [Tf{T^u) - V,E^{T^u))v{l) 

£eAMM\Ai 

:= Pi + P 2 + P 3 

with the interface region A^ := {£ G A, Rqm — Rbuf < Kl < Rqm + Rbuf}- 
To estimate Ri, we have from the expressions (5.3) that for any £ G A*^^\A^, 


(5.12) 


(5.13) 


(5.14) 




AQMuahuf 


R7(r^u) - V^E^iT^u) 

- E ( 

f-pSABUF 

- E 

^GpeASUF 

< Ce-^^BUF^ 


(Duo(£ -p)+ DT^u(£ - p)) - (Ruo(^ - p) + - P 


V, 


AQMuABUF 


Ip 


[Duo{£) + DT^u{£)) - V^^{Duo{£) + DT^u 


(5.15) 


wit 


A.2 


1 some constant p, where Theorem 3.1 (i) is used for the last inequality. Then we have from Lemma 
4 n _J D_^ 4 


that when Rbuf > ^ log Rqm and Rbuf > ^ log log Rmm, 

Ri < Ce-i^^^^\\Dv\\^2. 


(5.16) 
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To estimate P 2 , we have 


< \Pf^^{T^u) - PiiT^u)\ + \Vi£{T^u) - Vi£^{T^u)\ 

< C(^e-’’^BUF+ e-'r\p\ . \e^£ _ p) + 0 T'^u{£ - p)\Y) 

for £ G n and 


\P^{T^u)-Ve£^{T^u)\ 

< \P^{T^u) - Pf^^{T^u)\ + \Pf^^{T^u) - Pi{T^u)\ + \Ve£{T^u) - Ve£^{T^u)\ 

< C(^e-^^BUF+ Y e-'^\p\ ■ \e{e - p) + DT^u{i - p)\Y) 

^-peSflBUpWnAMM 


for £ G A^ n A'^^. Let A^^ := {i G A, -Rqm — -Rbuf < |^| < Rqm + 2i?BUF}- If P is satisfied, then we 
have from Theorem 13.31 and Lemma IA.2I that 

P 2 < C ^ (^e-^-f^BUF + |£|-m(fc+i)^ J' Y e-'^\P\ ■ \vii + p 
£eAi' IpIG^buf 


< C\\Dv\\i2 


logi^MM • .RbUF • (^Rqm ^ + -^QM • e ^^buf^ jf ^ 

• ( D“3A:-4/3 I n5/3 p-pilsuF^ 

-^BUF w riQjyj e j 


= 2 


if m = 3 


(5.17) 


If D is satisfied, then 


P 2 < C ^ |^e"^«BUF + |£|-(fc+i)^ Y • \v{e + p 

reA^' IpI<.Rbuf 

< Clogi^MM • 7 ?buf • + TZqm • ||L>r;||£2. (5.18) 

To estimate P 3 , let Pi{v) := Fi{Pxq + v) and £{v) := Yli^K {ReiR^o + v) — Ei{Pxq)). Define 

dTkj^i£{w) 


TkPi{w) = VeTk+i£{w) := 


dwf 


(5.19) 


Then, for any £ G A^'^\A\ we have 
\P^{T^u)-Vi£^{T^u)\ 

< \Pj^iT^u) - TkPeiuo + T^u)\ + |V,Tfc+iT(wo + T^n) - Vi£^iT^u) 

f Ce-P^BUFj£j-m jfp 

Ce-P^^v^\£\-‘^\og\£\ ifD ’ 


< 


where the same arguments as those in Lemma |B.3| are used to derive the last inequality. Then we 


have from Lemma 


A .2 


that when ifBUF > ^ log7?QM and T^buf > | loglogii’MM, 

P 3 < Ce" 4 ^BUFpy||^ 2 . 


(5.20) 
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Case P, m = 2 Case P, m = 3 Case D 

kp 

kp 

2 3 4 

1 2 3 

2 3 4 

1 2 3 

2 3 4 

12 3 

Rum 

error 

E-error 

p5 p7 

-^QM -^QM -f^QM 

p—3 p~5 

-^QM -^QM -^QM 

p—4 p~6 p~8 

-^QM -^QM 

p3 p5 p7 

■^QM -f^QM -f^QM 

p-4.5 p-7.5 p-10.5 

-^QM -^QM -^QM 

p—6 p~9 

-^QM -^QM 

Rqm ^qm -^qm 

p—1 p~2 p~3 

-^QM -^QM -^QM 
p—2 p~3 p~4 

-^QM -^QM -^QM 


Table 1: Choice of i?MM and error with respect to -Rqm for QM/MM schemes, with MM potential 
order k = for the energy based scheme and k = k-p for the force-based scheme. The energy error 
applies only for energy-mixing schemes. 


Taking (5.13), (5.14), (5.16), (5.17), (5.18), (5.20) and the estimates (4.25), (4.26) with order A;-|-l 


into accounts, we have the consistency 


\{T^{T^u),v)\<C\\Dv\U 




-|- e 4 


^BUF j£ p 


RJL log i?MM + RJm fog -^MM + e 4 -f^BUF if P) 


(5.21) 


when Rbvf > ^ logRqu and Rbvf > ^ loglogi^MM- 

4- Application of inverse function theorem: With the stability (5.12) and consistency (5.21), we 
can apply the inverse function theorem |23l Lemma 2.2] to obtain the existence of and the estimate 


n„ , + it P 

Du^ - T»u||^2 < <^ ) ( 

^ (^QM log -^mm + Rmm log Rmm + e"''-^BUF j if D 

with some constant k. This completes the proof. 


(5.22) 


□ 


6 Concluding remarks 

In this paper, we construct new QM/MM coupling algorithms for crystalline solids with embedded de¬ 
fects, based on either energy-mixing or force-mixing formulations. Unlike in commonly used QM/MM 
schemes, our approach does not employ “off-the-shelf’ interatomic potentials (or forces), but con¬ 
structs a potential (or force) specifically for the coupling with the QM model. The accuracy of our 
algorithms (with respect to increasing QM region size) is quantified by rigorous convergence rates. 

In the energy-based QM/MM coupling methods, with a given size Rqm of the QM region, we 
observe from Theorem 4.1 that one should take i?MM ~ -^qm (®-g-) fo ^1^® P) ^ = 2, ~ -^qm) 

and iiBUF ~ log Rqm to balance the errors. With these choices, we obtain the errors in Tablewritten 
in terms of 7 ?qM) dropping logarithmic contributions, and writing the order of expansion as k = kp. 

In our force-mixing QM/MM scheme, we obtain precisely the same rates and hence the same 
balance of approximation parameters, except that the order of expansion in the force is one less than 
that of the energy in our energy-mixing scheme. The rates are also shown in Table with k = kp. 

We note in particular that, for point defects, the QM/MM hybrid scheme achieves dramatic rates of 
convergence, already for a second order expansion of the site energies, respectively first order expansion 
of the forces {kp = 2,kp = I). By contrast, for dislocations, the second order expansion is no better 
than pure QM “clamped boundary condition” calculations (see [6l §4.2] and [9]). Only higher order 
expansions {kp > 3,kp > 2) of the site energy will give improved rates of convergence for hybrid 
QM/MM simulation of dislocations. 
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QM Region Radius [-RqmI 


QM Region Radius [-RqmI 


ATM 

QM/MM-Frc 
-k-k QM/MM-En 


Figure 6.3: Numerical verification of the convergence rates predicted in Theorems 4.1 and 5.1 (ATM 
denotes a pure QM scheme as described in [6]). The results are consistent with the theory, but the 
numerical rate for the energy is better than our analytical prediction. (See [9] for a similar gap in the 
theory). The inconsistency in the rates in the last data point in the energy error, and to some extent 
also visible in the displacement error for QM-MM-En, is likely due to a buffer radius that is chosen 
slightly too small for this level of accuracy. 


To limit the scope of the present work we will address the challenges in the implementation of both 
schemes in a separate article in full detail, but we present a preliminary numerical test. Using the TB 
toy model from m Sec. 5], the same simulation setup (2D triangular lattice with a di-vacancy defect), 
/ce = 2, fcp = 1, buffer radii Rb\jf = 1 + 0.61og(i?QM) and MM domain radii i?MM = 5 -Rqm + 2i?BUF 
we obtain numerical the results displayed in Figure |6.3[ This test should only be considered as a 
motivation for further study, but its implementation allows us to make the following observations: 

(1) A particular challenge in our schemes is the computational cost of higher order expansions, 
which is of the order ©((i^BUF)^™)• For example, taking only up to third neighbours in an FCC 
lattice (i?BUF/.RNN ~ 1-7, where i?NN is the nearest-neighbour distance) results in 42 neighbouring 
atoms, which would result in over 2M expansion coefficients at third order, and over 250M expansion 
coefficients at fourth order. We will exploit lattice symmetries to reduce the number of expansion 
coefficients that need to be calculated. The fact that the order of expansion is lower in force-based 
schemes, without loss of accuracy, is a significant advantage. 

(2) The computation of the fc-th order expansion of the site energies requires fe-th order pertur¬ 
bation theory (or, finite-differences). By contrast, the computation of forces and their derivatives can 
take advantage of the “2n-|- 1-Theorem”, hence expanding the forces is computationally much cheaper 
than expanding energies, even at the same order of expansion. An analogous comment applies to the 
computation of the QM region contribution to the hybrid forces or gradient of the hybrid energy. 

We conclude by commenting that, in view of the computational cost associated with Taylor expan¬ 
sions as site energies, alternative approaches may be required. Our analysis, or variations thereof, can 
then still be applied as long as the MM model is tuned to interact “correctly” with the QM model. 


Appendix A Interpolation of lattice functions 

For each n : A —>■ we denote its continuous and piecewise affine interpolant with respect to 7a by 
lu, and its piecewise constant gradient by VIu. We have the following lemma from [241125| . 

Lemma A.l. If v £ , then there exist constants c and C such that 

c\\Vvh2 < \\Dv\U2<C\\Vvh2. (A.l) 
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The following auxiliary results are useful in our analysis in that they sometimes allow us to avoid 
stress-strain (“weak”) representations of residual forces that we need to estimate. 

Lemma A.2. (i) If m = 2, then there exists C > 0 such that 

\v{I) — v{m)\ < C\\Dv\\( 2{1+ \og\l — m\) Vu G £,mGA. 

(ii) If m = 3, then there exists C > 0 such that, for each v G there exists Voa G such that 


V - Voahe < C\\Dv\\i 2 . 


Proof. The result is a straightforward generalisation of |24[ Proposition 12 (ii, iii)]. □ 

Lemma A.3. Let m = 2, 0<L<R and u : A —)• M satisfy v{£) = 0 V |£| > i?. // / : A —)> M satisfies 
|/(•^)| < then there exists a constant C such that 

Y. m • v{i) < Clog3/2 . \\Dv\U 2 V u G irH. (A.2) 

Proof. For simplicity of notation let r := |r|,r = r/r and v = Iv. Let R' > R, minimal, such that 
u = 0 in Bf^,. For each T G 7a, T C \ Bl we have 


< max \ f{i)\ Y^ b(^)l ^ f dr < (7 f r ^|Iu(r)| dr, 

c'eAnr ^sAnr ^gAnr ieAnr Jt 

Therefore, it follows that 


E /(^)-^(^) 

L<\l\<R 


< c 


|/u(r)| 




dr. 


(A.3) 


We have from the estimate 


that 


|u(r)| 


u r = 


r-R' 


d 


fR' 

< / |Vu(tf)| dt. 


dr = 


I- i-R' r-R' 

r“^|u(rf)| dr df < / / r~^ / |Vr(tr)| dt dr dr 

J Sm -1 J r=L Jt=r 




/S^-1 Jr=L 

r pB' ft 

|Vr(tr)| / r~^ dr dt df 

Jr=L 


s™-i Jt=L 

r r^' f 

i / t^/^|Vr(tf)|t“^/^ log — dt df 

s™-i Jt=L A 


^log^(i/A) dt df^ 

<Clog3/2 ||Vr||i2. 


Applying Lemma A.l, (A.3) and noting that log(i?Y72) < C completes the proof. 


□ 
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Appendix B Estimates of buffer truncations 


The following lemma can be proved by a calculation with the contour integration of the resolvent (we 
refer the detailed proofs to the arguments in [6l (3.10)-(3.12)]). 

Lemma B.l. Let i G M C N C A. Let y G Vm(A), 0 < j < n — 1, and p = (pi, • • • ,pj) G (M — £)T 
Then there exist constants C and rj, depending on m, such that 

\V^^p{Dy{L)) - V^^p{Dy{L)) \ < (B.l) 

with dist(l', n) = minfcgo{|^ “ Moreover, we have 

\V^^p{Dy{t)) - Vl^p{Dy{^)) \ < (B.2) 

A direct consequence of Lemma |B.1| is 

\Ve,p{Dy{i)) - Ve^^^{Dy{i)) \ < (B. 3 ) 

for y G Vm(A) and p = (pi, • • • , pj) G (A — ly with 0 < j < n — 1, max* \pi\ < i?BUF- 

For y G Vm, ^ G A^'^ U A^^ and i? > 0, we define V/^[Dy{i)) := [Dy{£)). Therefore, 

we have 


V/^p{Dy{i)) := 


if IpI > i? or £ + p G a n 


(^Dy{£)), otherwise 


The difference between Vg and can be estimated using Lemma 


B.l 


and p G (A — £y with max* \pi\ < R, then there exist constants C and p such that 

\Ve,p{Dy{£)) - Vl^p{Dy{i))\ < '^^O. 

The next lemma establishes the homogeneity of the site energy VP'. 


(B.4) 

If ii < 1^1 —i^DEF, 0 < j < n —1, 

(B.5) 


Lemma B.2. Let £,k G A and R < min{ 
for all IpI < R, then 


, \k\} - Rdef- If y,y' e Vo(A) satisfy Dpy{£) = Dpy'{k) 

(B. 6 ) 


Vg^{Dy{£)) = Vy^{Dy'{k)) and Vgl{Dy{£)) = V^^p{Dy'{k)) 
for p G (AZ™ n Br — Oy, max, |pj| < R. 


Proof. Using the condition {Dpy{£)}\p\^ji = {T)pp'(A:)}|p|<j:j and Theorem 3.1 (ii) with g{x) = x — 


x{£) + x{k), (iii) with G{x) = x — £ + k, we can derive Vf'[Dy{£)) = Vy^[Dy'{k)). Then the second 
part of ( B. 6 | ) is a direct consequence. □ 

Before the proof of (4.21), we need the following estimate for on the predictor, where 

the auxiliary site potential is defined by 


^/(g) := 


B^U-^DEF I 


g) 


IP 


K" (g) 


if \£\ < 3Rqm 
if \£\ > 3Rqm 


for g G (M"*) 




B.l 


Lemma B.3. Let Rb\jf > niax{^ logi?QM, ^ loglogi?MM}; where rj is the constant from Lemma 
If the assumption P or D is satisfied, then there exists a constant C such that 

Y, (<5U|UF(e(£)) -<5U/(e(£)),Du(£)) < Ce-^||Zlu ||,2 Vu G W^. (B.7) 

.feAMMuAFF 
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Proof. The left-hand side of (B.7) can be written in the form 


E E - Ve^p{Duo{£))) {v{£ + p) - v{i)) 

^g^MMuAFF p£A-e 

■= E (/^(^)-/™(^))^'(^)+ E (/^(^)-/''^^(^))^(^)> (B-8) 

^2QM-fiDEF<K|<3-RQM |^|>3 -Rqm 

where f^{i) and f^{i) are given in terms of V/^p{Duo{i)), and are given in terms 

of {Duo{£)). The precise forms of f^{£) and are not important, we can obtain from 

and the fact that they are given in terms of VPp{Duo{i)) and {Duo{£)) that 


Lemma 


B.l 


f^{e) - for _ i^DEF < \i\ < SRqm- 

When \i\ > 3-Rqm, we have 

^BUF(^) ^ 


/^(^) = 


E {vpypl{Duo{£ - P)) - Vpp^{Duo{l))) 

peAZ‘i-0, |p|<-Rbuf 

E {vPp^PDnoii - p)) - VPp{Duo{£))) 

p£AZ<i-0, |p|<4 


(B.9) 


and (B.IO) 


E 


p&AZ^-O, |p|<hl 


w 


Ve%ADuoii - p)) - Vp_pp{Duo{£ - p)) 


+ 


E 


peAi‘i-0, |p|<4 


w 


w 


v,AjDMe-p))-v,:i{DMi)) 


Lemma 


B.l 


:= fa{£) + Mi). 

and the definition of VP imply that 

fa<Ce-y\^\. 


(B.ll) 


(B.12) 


If P is satisfied, then we have from uo{i — p) = uq{£) = 0 and Lemma B.2 that 

h{l)=Q and = (B.13) 

If D is satisfied, we first consider the left halfspace ii < xi, in which case we can replace D by D. 
Let e = Duo{i) and suppressing the argument (i), we have from Lemma B.2 and the expansion 

I \A 

Vp{e) = Vp{0) + {5Vp{0),e) + - {l-t){6^Vp{te)e,e) dt with F = P/ or V|UF 

^ Jo 


that 


p,i 


Mi 


Mf) - /'*""(«) = E Ajo) - F#T(“) o-f'e 


+ lf(l-t)D-J[s^vfpte)-{S^V$^/{te))e,e\ dt := F, + F,. (B.15) 
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Using Lemma 3.1 and Lemma B.l we have 


which implies \fb{i) — /^^^(■^)| < Ce“^^BUF|£|-2 fg], < xi- For the right halfspace i\ > xi, we can 
repeat the foregoing argument to deduce 

\S {Mi) - \ < Ce-’?^BUF|£|- 2 _ 

Note that S is an 0(1) shift, which implies \fb{i) — < C'e“’''^BUF|£|-2 fgj. > fi. 

Taking (B.IO), (B.ll), (B.12), (B.13) and the above estimates for \Mi) ~ iiito accounts, 

we have 


< C + e-^^BUF|£|-2^ for \i\ > 3i?QM. 

Combining ( |B.8[ ), ( |B.9| ), ( |B.16[ ) and v G yields 


(B.16) 


with |f(.^)| < 


^SAMMuAFF RqM —^DEF<kl<.RMM 

Oe-^'^BUF if R(^y^ — i?DEp < |£| < 3i?QM 

C +e-’'^BUF|£|-2^ if |£| > 


Since i^suF > |logi2QM, we can write |f(l')| < Ce 2 ^buf|£| 2 ^ which together with Lemma 
completes the proof of (B.7) as i?BUF > ^ loglogi^MM- 


A.3 


□ 


Proof of (4.21). Denote g{£) = DT^rt(^) and suppressing the argument {£) in g{£) and e{£), we have 
from Lemma iB.ll and Lemma lB.31 that 

^ {6V^^^{e + g) - 5V/^{e + g),Dv{i)) = E (5F|UF(e) - 5U/(e), Du(£)) 

£g^MMuAFF ^g^MMuAFF 


+ E 


^gAMMuAFF ' 


\l-t) ((<5V|UF(e + tg) - 5V/(e + tg))g, Dr;(£)) dt 

< O (e- 2 ^BUF (b.19) 


Using (B.3), (B.5), ( B.19[ ) and Theorem |3.3[ we have 

Ti = E + -5U,(e + 5),Du(£)) 

^eAQM 

+ E (<^^|^^(e + <7)-5U/(e + 5),Du(£)) 

^gAMMuAFF 

+ E {SVt{e + g)-6V,{e + g),Dvi£)) 

£gAMMuAFF 

0(i?Q)Je"^'^BUF + g-I^BUF + 


< 


(B.20) 


which completes the proof since 7?Q^e ’?'^buf g gg^j^ fjg omitted compared with the middle 

term when T^buf > ^ logi^QM- 


□ 
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Proof of (4.31). Using T^u G AdmJ^ and the expansion 


^(e + g) = U(e) + (5U(e),5r) + - / {l-t){S‘^V{e + tg)g,g)dtwithV = ViOTVi 


BUF 


(B.21) 


for Rqm < I^I < Rmm + Rbvf, we have 

|5i| < {Ve{g + e)-Ve{e)-VP^^{g + e) + VP^^{e)) 

reAQM 

+ (sVe{e)-5VP^^{e),g) 


.Rqm< Kl <KmM+.RbUF 


1 

+ 2 


E 


(1 - t) {{6^Ve{e + tg) - S^VP^^{e + tg))g,g) dt 


.RQM<h|<J?MM+JJBUF 

:= Sf + S'l + SI. 


(B.22) 


The first and third group of (B.22) can be estimated by Lemma B.l 


5 ?| < Ci?QjJe"^^BUF and \S1\ < Ce-'^^^^^\\g\\%. 


(B.23) 


Using Lemma B.3, we can obtain the estimate for the second group: 

\Si\ < Ce^^^^uF when i^suF > - log log i^MM, 

V 


which together with (B.22) and (B.23) completes the proof of (4.31). 


□ 


Appendix C Stability of force-mixing methods 


Here, we establish the result that the energy-mixing Hessian and force-mixing Jacobian are “close”. 
This result is reminiscent of similar results in the context of atomistic/continuum blending [19], but 
the proofs are not closely related. 


Proof of (5.11). Let := 
Denoting := 


{£ G A, Rqm — .Rbuf < l-^l < Rqm + .Rbuf} be the interface region, 
and := {E^'"^{Pxo + v) - Ef^^{Pxo)) with = 
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we can split 


{{6T^{T^u) - 6^S^{T^u))v,v) = '^{6T^{T^u) - 6VeS^{T^u),v)v{e) 


£eA 


{6Tf{0)-6V,£'^{0),v)vii) 


£eAZ" 


+ Y 1 {SJ^fiT^u)-5ViS^(T^u),v)v{e) 

£eAQM\Ai 

Z (5-F;’(0)-<5V,^(0),n>n(^) 

^e^^'^nSKQM-flBUF 

+ ^ ({(ijf (T^n) - 6TX0),v)v{£) - {SVeS^{T^u) - 6Ve^{0),v)v 
£eAi 

+ ^ (^{6T^{T^u)-6T^{0),v)v{£)-{6Ve£^{T^u)-6Ve^{0),v)v 

£eAMM\Ai 

:= Qi + Q 2 + Qs + Qa: + Qs- 

Estimate for Qi: Using Theorem [3.2| (ii), (iii) and m Lemma 3.4], we can rewrite Qi as 


^1= E E Dpviif{A^-A^)Dpvii) 

p&AZ'^-O 


(C.l) 


(C.2) 


with Ap,A^ G Ap = —2-^o,p(0) A^ = —5^^op(d)- This crucially uses that fact that we 

apply the force approximation at every lattice site, which makes it conservative. Note that Lemma 


B.l and Theorem 3.1 (i) imply 


1^-^ _ < (jQ-n{RBV¥+\p\) 


which together with (C.2) yields 


IQil <Ce-^^^'^^\\Dv\\l. 

7 

Estimates for Q 2 , Q 3 : Lemma A.2| and v G imply that 

||n||£oo < CllZlnll^^ (1 + log I^mm) ) if m = 2, and 
\\v\\e 6 < C\\Dv\\i 2 , ifm = 3. 


(C.3) 


(C.4) 


IQ 2 I < C e-^^^^^\Dvii)\^\v{i)\ 


< C 


Using (C.4), Lemma B.l and Theorem |3.1| (i), we have 

3-piJBUF I 

teAQM\Ai 

-Rqm • logi?MM • e“^-^BUF||/)^||2^ if m = 2 

< Ci^QM • logi?MM • . 

Analogously, we have 

m < CRqm • logRuM • e-^^^^^\\Dv\\%. 


if m = 3 


(C.5) 

(C.6) 
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Estimate for We can rewrite Qi as 

Qi = Y^ - 5Ff^^{T^u),v)v{l) + ^ {5Ff^^{T^u) - dF^uo + T^u),v)v{l) 

reAi 

- {5Vi£^{T^u) - 6Vi£^^^{T^u),v)v{l) - Y _ sVi^{uo + T^u),v)v{l) 

i&F reA^ 

+ ^ {6F^{uo + T^u) - SF'fiQ), v)v{£) - Y (<^Vr^(uo + T^u) - 6Vi^{0),v)v{i). 

£eAi £eAi 


Using and ^^1 < Ce“^(l^®uF-|p||+|p|+l?l+ICl) 


U-p,pCC 


to Lemma B.l), the last line an be bounded by 


(this is proved analogously 


E E ( -p) + -p))- (0)) 

teAipe^-Ab |p|<-Rbuf 

- {ivf4F'^‘’\e{e -p) + DT“s((> - p)) - , 0v{t - p)) v{l) 


= E E /‘((<5^**7-''>{t(e(«-p) + Dr«ii(f-rt)) 

l&A^ |p|<J?BUF 

- -p'^ + ^T^u{i - p)))) (e(£ -p) + DT^u{l - p)),Dv{l - p)) dt v{i) 


< C E • |e(^) + • |Du(£)|^ • ( Y • |u(£ + p)|). 

£eAi' IpI<-Rbuf 


with := {I G A, max{0, i?QM — 2 -Rbuf} < Kl < -Rqm + 2-Rbuf}- Using Theorem 3.3 and (C.4), 
we have 

\Qi\ < C (|e(£) + DrHu(£)|^ + e-^^BUF^ . ^ g-AH . |^(£ + p)|^ 

£eAi' IpIG^buf 


— y'^BVF ■ < 


log i?MM (.Rqm ^ »?^BUF^ if P with m = 2 
+ e“’'^suF if p with m = 3 

logRMM{Rfil~^^^ + if D 

< C<u'f • log ^mm • (^QM+ e-^^BOF) 


\\- 


(C.7) 


Estimate for Q^: Let Fi{v) := F^{Px + v) and £{v) := Yle&A {Ri{R^o + v) — Ei{Pxo)), then 
Fi{v) = Vi£{v). Further, we define 


fkPeiw) = Vefk+i£{w) := where (C.8) 

OWi 

fk+i£{w) = Tk+i£{w) - T(0) - (0), w). 
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Then, for any i G we have 


{8Ff{T^u) - (^Jl’(O)) - {5Vi8^{T^u) - <5V^^(0)),n 


< 


<5J7(T^n) - <57^(0) - 5TkTi{uo + T^u),v 


+ 


{ SVeS^iT^u) - SVeS^iO) - 6ViTk+iS{uo + T^u), 


< (7g“^^BUF E • \Dv{i-p)\y 

\p\<Rbvf 


\£- p\-^ 

\i- p\~^ log\i- p\ 


if P 
if D 


Summing over £ G A'^^ \ A\ and applying (C.4) we obtain 

IQsl < C log2 ^ • e-^^^^^\\Dv\\%. 

Aqm ^ 


(C.9) 


(For case P one obtains log instead of log^ qualitatively the same as replacing 

the unknown exponent p with p/2, hence we ignore this small improvement.) 


Combining ( |C.1| ), ( |C.3 ), (C.5), ( C.6| ), ( |C.7 ) and (C.9), we finally deduce that 




provided that Rqm > log'^7?MM, 7 ?buf > |logi?QM, 7 ?buf > |loglogi?MM- This completes the 
proof. □ 
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